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FROZEN GAUSSIAN APPROXIMATION FOR HIGH FREQUENGY WAVE 
PROPAGATION IN PERIODIC MEDIA 

RICARDO DELGADILLO, JIANFENG LU, AND XU YANG 


Abstract. Propagation of high-frequency wave in periodic media is a challenging problem due to the 
existence of multiscale characterized by short wavelength, small lattice constant and large physical 
domain size. Conventional computational methods lead to extremely expensive costs, especially in high 
dimensions. In this paper, based on Bloch decomposition and asymptotic analysis in the phase space, 
we derive the frozen Gaussian approximation for high-frequency wave propagation in periodic media 
and establish its converge to the true solution. The formulation leads to efficient numerical algorithms, 
which are presented in a companion paper [ 3 . 


1. Introduction 

We are interested in studying high-frequency wave propagation in periodic media. A typical example is 
given by the following Schrodinger equation in the semiclassical regime with a superposition of a (highly 
oscillatory) microscopic periodic potential and a macroscopic smooth potential, 

( 1 . 1 ) + V{x/+ U, a; G 

where V and U are smooth potential functions, V is periodic with respect to the lattice V{x + 6^) = 
V(x) for any a: G and {e^, i = 1,2, •• • ,d} is the standard basis of Here e <C 1 is the rescaled 
Planck constant, is the wave function, and d is the spatial dimensionality. 

The equation (HU can be viewed as a model for electron dynamics in a crystal, where V is the effective 
periodic potential induced by the crystal, and U is some external macroscopic potential. Notice that we 
have identified the period of V(x/e) and the “semiclassical parameter” in front of the derivative terms. 
This parameter choice gives the most interesting case as e ^ 0 [5]. 

The mathematical analysis of this work is motivated by the challenge of numerical simulation of (HU 
when e is small. In this semiclassical regime, the wave function becomes oscillatory with wave length 
0{e). This means a computational domain of order 1 size contains 0(1/e) wavelengths, and each of 
them needs to be resolved if conventional numerical methods are applied. For example, even for the 
simplest case V = 0 (no lattice potential), a mesh size of 0(e) is required when using the time-splitting 
spectral method H] to compute (HU directly; an even worse mesh size of 0(e) is needed if one uses 
the Crank-Nicolson schemes [26] or the Dufort-Frankel scheme [27|. Besides, the presence of non-zero 
lattice potential introduces further difficulties which restrict the mesh size to be o(e) in the standard 
time-splitting spectral method [T] ■ Special techniques using Bloch decomposition are needed to relax the 
mesh size to be of 0(e) [12UI4] . Moreover, in these methods, a large domain is demanded in order to 
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avoid the boundary effects. Therefore the total number of grid points is huge, which usually leads to 
unaffordable computational cost, especially in high (d > 1) dimensions. 

An alternative efficient approach is to solve cu asymptotically by the Bloch decomposition and 
modified WKB methods Em 13, which lead to eikonal and transport equations in the semi-classical 
regime. An advantage of this method is that the computational cost is independent of e. However, the 
eikonal equation can develop singularities which make the method break down at caustics. The Gaussian 
beam method (GBM) El] then introduced by Popov to overcome this drawback at caustics. The idea 
is to allow the phase function to be complex and choose the imaginary part properly so that the solution 
has a Gaussian profile; see [HHinilMlEIlEE] for recent developments. Similar ideas can be also found 
in the Hagedorn wave packet method ElEI- Unlike the geometric optics based method, the Gaussian 
beam method allows for accurate computation of wave function around caustics [611^. But the problem 
is that the constructed beam must stay near the geometric rays to maintain accuracy. This becomes a 
drawback when the solution spreads 

The Herman-Kluk propagator umim was proposed for Schrodinger equation without the oscillatory 
periodic background potential. The method was rigorously analyzed in [351136] and further extended as the 
frozen Gaussian approximation (FG A) for general high frequency wave propagation in [23H25j . The FGA 
method uses Gaussian functions with fixed widths, instead of using those that might spread over time, 
to approximate the wave solution. Despite its superficial similarity with the Gaussian beam method, 
it is different at a fundamental level. FGA is based on phase plane analysis, while GBM is based on 
the asymptotic solution to a wave equation with Gaussian initial data. In FGA, the solution to the 
wave equation is approximated by a superposition of Gaussian functions living in phase space, and each 
function is not necessarily an asymptotic solution, while GBM uses Gaussian functions (called beams) in 
physical space, with each individual beam being an asymptotic solution to the wave equation. The main 
advantage of FGA over GBM is that the problem of beam spreading no longer exists. 

In this paper, we extend FGA for computation of high-frequency wave propagation in periodic media. 
We mainly focus on the derivation of an integral representation formula of FGA in the phase space and 
establish the rigorous convergence results for FGA. While the FGA works for general strictly hyperbolic 
equations, we focus in this paper the case of semiclassical Schrodinger equation with periodic media (EU. 
The computational algorithm and numerical results will be presented in a separate paper E] • The rest of 
the paper is organized as follows. We first recall the Bloch decomposition of periodic media and introduce 
the windowed Bloch transform in Section [2] In Section El we present the formulation of FGA for periodic 
media and the main convergence result. The proof of the main result is given in Sectional 

Notations. The absolute value. Euclidean distance, vector norm, induced matrix norm, and sum of 
components of a multi-index will all be denoted by | • |. We will use the standard notations S, C°°, and 
for Schwartz class functions, smooth functions, and compactly supported smooth functions, respectively. 
We will sometimes use subscripts to specify the dependence of a constant on the parameters, for instance, 
notations like Ct to specify the dependence of a constant on a parameter T. 


2. Bloch decomposition and windowed Bloch transform 

The frozen Gaussian approximation for periodic media relies crucially on the Bloch decomposition to 
capture the fine scale (0{£) spatial scale) oscillation. First we briefly recall the well-known Bloch-Floquet 
decomposition for Schrodinger operators with a periodic potential. 
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Consider a Schrodinger operator 

(2.1) if = -iA +1/(01), 

where the potential V is periodic with respect to the lattice Z'^. We denote F the unit cell of the 
lattice: F = [0,1)'^. The unit cell of the reciprocal lattice (known as the first Brillouin zone) is given by 
F* = [—7r,7r)‘^. It is standard (e.g., M) that the spectrum of H is given by energy bands 

CXD 

spec(ii) = u u Enii), 

where for each ^ G F*, {En(^)} are the collection of eigenvalues (in ascending order) of the operator 

with periodic boundary condition on F. The Bloch waves are the associated eigenfunctions: For each 
band n and ^ G F*, it solves 

(2.2) ff^UnU,-)=Er,U)UnU,-)- 

with periodic boundary condition on F, where £ serves as a parameter in the above equation. Un(^, •) is 
normalized that 

(2.3) J dx = 1. 

We extend Un(^,x) periodically with respect to the second variable, so it is defined on F* x R'^. We will 
also write = Un(^, ■) when the former is more convenient. 

These Bloch waves generalize the Fourier modes (complex exponentials) to periodic media (see for 
example discussions in m)- In particular, for any function / G L^(R‘^), we have the Bloch decomposition 

(2.4) f{x) = Y. I, =*=)e‘« "(e/)„(0 de 

where the Bloch transform B : T^(R'^) —>■ L^(F*)^ is given by 

(2.5) (g/)n(|) = , f Mn(l,2/)e”'^'"/(y)dy. 

(271)'*/^ jRd 

As an analog of the Parseval’s identity, we have 

p OO p 

(2.6) / |/(a;)|' da; = ^ / \{BfU^)f 

As suggested by (12.41) and (12.51) . we introduce the notation 11 to denote the phase space corresponding to 
one band (F* is viewed as a torus, i.e., periodic boundary condition is assumed on F*) 

(2.7) n:=W^ xT* = {{x,^) \xgW^,^G F*}. 

Correspondingly, we will use the notation (q, p) for a point in fl. 

For later usage, we define the Berry phase An for the Bloch waves, 

(2.8) A„(|) = (Mn(l,-),iV^'Wn(l,-))i2(r)- 

The normalization condition (12.31) implies /I„(^) is always a real number. We should be cautious about 
one subtlety though as the eigenvalue equation (12.21) and the normalization only dehne Un(^, •) up to a 
unit complex number, in particular, for any function ip periodic in F*, 

iv(C) 


(2.9) 


Vn{^,x) = e 


Un{^,x), {x,^)Gn, 
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also provides a set of Bloch waves. This is known as the gauge choice for the Bloch waves. However, 
different gauge choice gives different values of An{^) and even causes trouble if ip is discontinuous. While 
for the analysis, it suffices to assume smooth dependence of on ^ (which is possible as the n-th band 
is separated from the rest of the spectrum), this gauge freedom makes numerical computation nontrivial. 
We will further address this by designing a gauge-invariant algorithm in a companion paper on the 
numerical algorithms. 

Differentiating (12.211 with respect to ^ produces 



Taking inner product with Un{^-, •) yields 


( 2 . 11 ) 




Differentiate (|2.10l) with respect to ^ again gives 


( 2 . 12 ) x) + 2 (-iVa; -f |)V^u„(|, x) + u„(|, x)I 


^^™(l)v|u„(|, x) + 2 V^£'„(|)V^u„(|, x) + En{i)Vlun{$, x). 


Taking inner product with Un{£,, ■), one gets 


(2.13) {Unii, ■), -iVa;V^M„(|, a;)) -k |('u„(|, •), •)) + 



These identities (12.111) and p.l3l) will be useful later. 

We shall now introduce the windowed Bloch transform. This is an analog of the windowed Fourier 


transform (also known as the short time Fourier transform) widely used in time-frequency signal analysis. 


Definition 2.1. The windowed Bloch transform W : L^(IR'^) ^ T^(D)^ is defined as 



(2.14) 


where Gq.p is a Gaussian centered at (q,p) G D, given by 



(2.15) 


The adjoint operator W* : L^(D)^ —>■ T^(R'^) is then 



n—1 


Proposition 2.2. The windowed Bloch transform and its adjoint satisfies 


>V*W = ldi2(„.). 


(2.17) 


Remark. Similar to the windowed Fourier transform, the representation given by the windowed Bloch 


transform is redundant, so that WW* 7 ^ ldi 2 (Q)N. The normalization constant in the definition of W is 
also due to this redundancy. 
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Proof. Fix a / S L^(]R“), by definition, we have 

2^dj2 


^dj Z ^ r r 

{W*Wf){x) = Xj JJ^^n{p,x)Gq^p{x){Gq^pUn{p,-), f) dq dp 


c^d/ Z ^ f f f — 

(27i-)3ci/2 Gq.p{x)Gq.p{y)un{p, v) f {v) dy dq dp. 


Let us integrate in q first. 


Gq^p{x^Gq p(^y^ dq — c 


— Jp-{x-y) / ^-\x-q\ /2-\y-q\^/2 


dq 


= eip (^-j/)g-N-y|V4 / exp|- 
JR-i 

_ .j^d/2^ip-(x-y)^-\x-y\^ /4^ 


q- 


x + y 


dq 


Hence, denoting f^{y) = e we have 

1 OO p p 

(W*yv/)(a;) = —— V / / M„(p,a;)e*P'(“="2/)e-|^-y|V4y^(p^y)j(y)dydp 

(27r)<i 


f f Un{p,x)pPd^ ^^w„(p,y)/^(y)dydp 


iMl 


= /.(a:) = /^f{x)=f{x). 


□ 


The previous proposition motivates us to consider the contribution of each band to the reconstruction 
formulae (I2.17L This gives to the operator Hj^ : L^(]R^) —?> L'^(W^) for each n G N 

nd/A r p 

(2-18) {U^f)ix) = JJ^Unip,x)Gq^p{x){'Wf)n{q,p)dqdp. 

It follows from (12.171) that = IdL 2 (Rd), while H^ is not projection due to the redundancy of 

windowed Bloch transform. 


3. Formulation and main results 


Let us start with fixing some more notations. We will switch between physical domain and phase space 
in the FGA formulation. For clarity, we will use x,y G as spatial variables, (g,p) G as phase 
space variables. The capital letters X and Y are shorthand notations for X = a:/e and Y = y/e. 

We define an effective (classical) Hamiltonian corresponding to each energy band by 


(3.1) 


hn{q,p) = En{p) + U{q). 


The associated Hamiltonian flow Kn{t) = q,p), Pn{t, q,p)) solves 


(3.2) 


dQ 

dt 

dt 


n 


n 


VPn {Qnj Pn)i 


-'^Qr.Qn^Pn) 


on H with initial conditions Q„(0, q,p) = q and P„(0, q,p) = p. 

From now on, we will use the short hand notation ((5„,P„) for ((5„(t, q,p), P„(t, g,p)). For the 
long time existence of the Hamiltonian flow (13.2L we will assume that the external potential U{x) is 
subquadratic as below. 
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Definition 3.1. A potential U is called subquadratic, if ||9“C/(a;)||^oo is finite for all multi-index |a| > 2. 
Remark. As a result, since the domain F* for p is bounded, the Hamiltonian is also subquadratic. 


The frozen Gaussian approximation will be formulated by the following Fourier integral operator. 

Definition 3.2. (Fourier Integral Operator) For u G x H, C) and ip G 5(M‘^,C) we define the 

Fourier Integral Operator with symbol u by the oscillatory integral 

(3.3) [I^{u)ip\{x) = \ [[ [ ei‘^^*’^’y’‘>^P^u{x,y,q,p)q:}{y)dydqdp 

[Zire) / JJaJRd 

where the complex valued phase function ^{t,x,y,q,p) is given by 

(3.4) <^{t,x,y,q,p) = S{t, q,p) - p ■ {y - q) + P ■ {x - Q) + ^\y - qf ^1® - Q|^ 
and S{t,q,p) is a real-valued action function associated to n satisfying 

(3.5) VgSit,q,p) = -p + VgQ-P, VpS{t,q,p) = VpQ-P. 


Note that if K{t) = /c„(t), the action Sn{t, q,p) can be obtained by solving the evolution equation 


(3.6) 


^=Pn- ypMQn, Pn) " K{Qn. Pn), 


with initial condition 5„(0,q,p) = 0. 

We are now ready to formulate the frozen Gaussian approximation. The FGA approximates the 
solution of the Schrodinger equation (ED on the n-th band to the leading order by 

(3.7) (^an,oit,q,p)uniPn,^)un{p,^)jtpo 

where "^q is the initial condition. More explicitly, at time t, i/’pqa is given by 

(3.8) #GA(i,®) = (27r£)3^/2 ix)u„{Pn, x/e) 

■ (Gq.pU„(p, •/£),'0o)dgdp. 

Here and in the sequel, we use the short-hand notation for Gaussians with semiclassical scaling 


(3.9) 


GqJx) ■■= exp - 


\x-q[ 

2e 


-hi 


■ P- [x-q) 


where the subscripts {q,p) indicate the center of the Gaussian in phase space. Note that the semiclassical 

K-P? , .9 • (I -P)' 


Fourier transform of G^p is 


(3.10) 


SiPS) = dx = exp (^- 


2e 


-hl- 


Gorrespondingly, the semiclassical windowed Bloch transform 

2C^/4 


i^(H)^ is defined as 


2Ci/4 f 

(3.11) (W^/)„(g,p) = (27r£)3^/4 ■/^)Gq,p, f) = (27r£)3rf/4 ^/^)Gq,p{x)f {x) dai. 

for each n G N with semiclassical scaling 


Similarly we also have the operator 

2^/4 


(3-12) {U^''^f){y) = ^ 2 jg) 3 d /4 y/s)Gl^^{y){yV^ f)u{x, 0 <ix d|. 


ft follows from (12.171) and a change of variable that J2n 
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The only term in (13.8|) that remains to be specified is the amplitude an,o{t, Q,p)- It solves the evolution 
equation 

(3.13) dtttn,0 = -ian,oAn{Pn) ■ VU{Q„) + ^an,oti (dzPnV^E{Pn){Zn) 

- ^a„,o tr . 

with initial conditions a„^o(0, q,p) = 2"^/^ for each {q,p) and we recall that = {un{^, •), iV^it„(^, •)) 

is the Berry phase. Here the matrix Z associated with the Hamiltonian flow K„(t) is defined by 

(3.14) Zn{t,q,p) := d^{Q^ + iPn), 
where := dq — idp. 

We now state the main results of this work. 


Theorem 3.1. A ssume that the n-th Bloch band if„(^) does not intersect any other Bloch bands for all 
^ G r* and the Hamiltonian hn{x,^) is subquadratic. Let be the propagator of the time-dependent 
Schrodinger equation CH) with initial condition € i^(R‘^). Then for any given T, 0 < t < T and 
sufficiently small e, we have 


(3.15) 


sup 

0<t<T 


-I‘'{a„pUn{Pn,x/e)un{p,y/e))% 


L2 


< Ct.u £ \m 


0||_L2- 


Remark. Note that the FGA solution approximates the time evolution of which is the n-th band 

contribution to the initial condition in the reconstruction formulae ()2.17p . In particular, if the initial 
condition is concentrated on the n-th band in the sense that the theorem states that the 

solution to dni) is approximated by the FGA solution with 0(e) error. 


Remark. We can also construct higher order approximations by replacing the term anpUn{Pn, x/e) with 
an e-expansion of the form bnp + + e^&n ,2 -I- ... -I- £^~^bn^N-i where bnp = an,oUn{Pn, x/e). This 

will give an approximate solution 1° 0{e^) accuracy. In this paper we shall focus on the first order 

approximation and omit the formulation and proof for higher orders. 


Remark. Let us also remark that while we take the more explicit approach of using Bloch waves in a 
modified FGA ansatz for periodic media, as in (13.81) . The same approximation can be also derived by 
first projecting the whole Schrodinger equation using a super-adiabatic projection as developed in [2^1^ 
and then apply the frozen Gaussian approximation to the resulting dynamics. We will not go into the 
details in this work. 


The proof of Theorem 13.11 is given in Section S) By linearity of (ll.l|) , we have the following more 
general statement, as an easy corollary from Theorem 13.II 

Theorem 3.2. Assume that the first N Bloch bands n = 1, - ■ ■ do not intersect and are 

separated from the other bands for all ^ G F*; and assume that the Hamiltonian hn{x,^) is subquadratic. 
Let be the propagator of the time-dependent Schrodinger equation (II.fl) with initial condition S 
Lf{W^). Then for any given T, 0 < t < T and sufficiently small e, we have 


(3.16) sup 

0<i<T 


N 


n —1 


L2 


{anfiUn{Pn,x/e)Un{p,y / £))lp[ 

< Ct,n ^ 11 ^ 011^2 


N 




L2 
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Proof. Taking the short-hand notation i/jg ^ and = I^(an,o'U'n{Pn,x/e)un{p,y/e)), we 

have 


N 


- E 


L 2 


N 


E 'i’ln - E nui’l 


\n—l 
/ N 


L 2 


N 


E ^o,n ]+^n E ^o.„ - E 


< 


^n=l 
^ N 


\n—N+l 


L 2 


iV 


E ^O.n - E 

\n=l / n=l 

^ ^ „ „ 

< 2^ C'T,ne 11 11^2- 


L 2 


E 


0 ,n 


^n=Ar+l 


L 2 


137151 


E « 

n=A^+l 


< Ct,NS 11^0 IIl2 + 


N 




n—1 


L 2 


L 2 


□ 


4. Analysis of frozen Gaussian approximation in periodic media 


4.1. Initial condition. Let us first study the initial condition for the frozen Gaussian approximation. 
At time t = 0, observe that by setting f = 0 in (13.8|) we have 

2^/2 


2^/2 PP 

i’FGA,= (27r£)3^/2 JJ^^ri{p,x/e)G‘^g^p{x){G‘^gpUn{p,-/£),i^Q) 


dqdp = U^’^ro 


by definition of the operator Hence, the FGA solution matches at t = 0. 


4.2. Estimates of the Hamiltonian flows. To control the error for t > 0, we collect here some 
preliminary results on the estimate of quantities associated with the Hamiltonian flows. We will assume 
throughout the rest of the paper that the assumptions of Theorem 13.II hold for a fixed Bloch band n. 
The following notation is useful in the proof. For u G C), we define for k G N, 

(4.1) Mk[u] = max sup \dq'’ dp^ u{q,p)\ 

l“gl + l“pl<fc (q,p)en 

where and ap are multi-indices corresponding to q and p, respectively. 


Definition 4.1. (Canonical Transformation) Let k : —>■ be a differentiable map K{q,p) = 

{Q{q,p), P{q,p)) and denote the Jacobian matrix as 

(. 2 ) ({dqQV{q,p) {dpQ)^{q,p)\ 

^ - ^ > \{dqPnq,p) idpP)^iq,p)) ■ 

We say k is a canonical transformation if F is symplectic for any {q,p) G R^'^, i.e. 


(4.3) 





It is easy to check by the definition that the map K„(t) : R^"^ R^'^ defined by (t/jp) —> {Q^{t,q,p), Pn{t,q,p)) 

solving (IX^ is a canonical transformation. 
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Proposition 4.2. We have for all k > 0 


(4.4) 


sup Mfc [Fn{t)] < oo sup Mk 

te[0,T] t6[0,T] 




< oo. 


Proof. Differentiating Fn{t,q,p) with respect to t gives 
(4.5) 


d , . / dpdqhn dpdphn \ p /, N 

—F„(t,q,p)= „ o a . \Fn{t,q,p). 

at y—dQOQhn —Oqdphn I 


By our assumption that U is subquadratic on and since En € C°°(r*), there exists a constant C 
independent of {q,p) such that 


(4.6) 


— \Fr,{t,q,p)\ = 


Fn{t,q,p)\ < C\Fn{t,q,p)\ 


dpdqhn dpdphn 
^—dqdQhn —dqdphn 

with |i^„(0)| = |Id 2 d|- By an application of Gronwall’s inequality, we obtain 

(4.7) \Fr,it)\<e^\^\. 

Differentiating (14.51) with respect to (q,p) yields 

(4.8) 43.°’«'A(i.9.p)= E (°f)(7y^'^’(vfk 

dt \PvJ \-dQdQhn -dqdphrk) 

xdy-^’’d;^-^^Fr,it,q,p). 

Our estimate now follows by induction. □ 

Recall that the matrix Zn{t,q,p) is defined by 

(4.9) Zn{t,q,p) := 9^ (Q„(f,q,p)+iP„(t,q,p)) = {dq - idp) {Q^{t,q,p) + iPn{t,q,p)). 

We have the following. It follows the same proof of [331 Proposition 3.5], which we reproduce here for 
completeness. 

Proposition 4.3. Zn{t,q,p) is invertible for {q,p) G D. Moreover, for each fc € N, 


(4.10) 


Mk 


{Zn{t)y 


< oo. 


Proof. Zn{t, q,p) inherits the property that Mk{Zn{t, q,p)) < oo from the same estimate for F„(t, q,p). 
Moreover, we have 


i{Zn)*{t,q,p) = (^ildd Idd) {Fnf{t,q,p) f 


Id, -ild,\ -ild. 

Id, Id, j Ud , 


(4.11) 


(ild. Id,) ((F„)^(F„))(t,9,p) 


+ (ild. Id,) {Fnf{t,q,p) 


-ild,\ 

. / 

0 -ild,^ 


IV 0 

-ild,^ 


(ild. Id,) ((F„)^F„)(t,g,p) ' ^ +2Id,. 


Id, 

This calculation shows that, since {FnftfyFn{t) is semi-positive definite, for any v G 

v*Z,,{t)iZ^{t))*v>2\v\7 


(4.12) 
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Therefore Zn{t,q,p) is invertible and det[Zn{t)) is uniformly bounded away from 0 for all q and p, so 
by representing {Zn)~'^{t,q,p) by minors, Mk{{Zn)~^{t,q,p)) < oo, as Mk{Zn{t,q,p)) is. □ 

Proposition 4.4. For each k €N, 

(4.13) sup Mfc [Un{Pn,x)] < OO. 

te[o,T] 

Proof. Un{Pn,x) is smooth on the compact set T* x T since the n-th band is separated from the rest of 
the spectrum (see e.g., [MJ Sec XIII.16]). Thus Un{Pn,x) is uniformly bounded on T* x T and hence 
r* X due to periodicity. We also see from Proposition 14.21 that the derivatives of u„(P„,a;) are also 
bounded. Thus, Mfc[u„(P„, a:)] < oo for any finite time t. □ 

4.3. Higher order asymptotic solution. To prove the theorem, we will need to construct a solution to 
the Schrodinger equation that is accurate up to 0{e‘^). The construction is based on matched asymptotic 
expansion. Let us fix a band n G N and consider the ansatz 

(4.14) V'FGA.oo = (G^_pU„(p,-/e),V'o) dqdp, 

where the coefficient b assumes the asymptotic expansion 

OO 

b%t,X,q,p) :='^e^bj{t,X,q,p) 

3=0 

(4.15) 

+ q,p)Un{Pn^ X) + Q,P)) 

oo 

+ G^(an, 2 (t, q,p)un{Pn-, X) + X, g, p)) + ^ e^bj(t, X, g, p) 

3=3 

To determine the terms in the expansion, we will make use of the following Lemma. 


Definition 4.5. For / = f{t,x,y,q, p) and g = g(t,x,y,q, p) such that for any t and x, 

fit, X, .,., •), 5(t, ai, •, •, •) G L^iR‘^-, SiR‘^ x T*)), 

we say that / and g are equivalent for the n-th Bloch band, denoted as / 9 if for any t > 0 and 

(4.16) [[ [ (/-g)(t,ai,y,g,p)G^ p(y)4'o(y)dydgdp = 0. 

J Jn Jr^ 

Lemma 4.6. For any d-vector function v{y,q,p) such that each component is in L°°(R‘^; <S(R‘^ x T*)) 

(4.17) v{y,q,p) ■ {x - Q„) -ed^ ■ (v^”^), 

and for any dx d matrix function Miy,q,p) such that each component is in L°“(R‘^; <S(R‘^ x T*)) 

tr (M(y, q,p){x - QJ'^) tr [d^Q^MZ-^) - etr [d^Mix - QJZ-'^ + Mix - Q^)d^Z-^) 

=Gtr id,Q^MZ-^)+eHT {d, {d,MiZ-^f) + 8, {Md,Z-^) Z-^) . 
Higher order terms can be obtained recursively. In general we have for any multi-index a that |a| > 3, 
(4.19) ix-QX -nO{e\-^\). 

Proof. The proof of lemma |T6] is essentially the same as in Lemma 3 of [36] and thus is omitted here. □ 
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We now substitute (14.141) into the Schrodinger equation. For this we first compute the time and space 
derivatives on iPfga oo ■ 


(4.20) ie^tV'FGA.cx. = (27r£)3^/2 J Xj - idtQJ ■ {x - QJ) 6^} x 

X •/e),'i/'o) dqdp. 


(4.21) Ai/ipQA^oo — 


(27re)3'i/2 


-2(-iVx + Pnfb^ - (Vxfo" + ib^Pn) ■ {x - Q„ 
+ ^b^\x-QJ^ -^eb^d 
X G^^_p„e*'®"/^(G^ pU„(p, •/e),'0o) dqdp. 


Hence, after rearranging terms, we arrive at 
(4.22) 

{iedt + - y(X) - G(a;))#GA.oo = 


1 


(2^£)3d/2 JJ^ 


-2(-iVx + Puf - V{X) - U{x) - dtSr, 


+ e{idtb^ - -¥d) - [(Vx6^ + d>^Pn) + {dtPu - - QJ+ 

+ \\^~ Qnl'^b^ + Pn ■ 5tQ„6®|G^^ (G® pU„(p, ■/e),'ipo) dqdp. 


Define 


f{t,x,y,q,p) = 


(4.23) 


- 2 (-iVx + Pnf - V{X) - U{x) - dtS„ 


+ e{idtb‘^ - -b^d) - [{Vxb^ + ib^Pn) + {dtPn - idtQnW] ' “ Qn)+ 

+ ^\x- . dtQ^bAunip, Y), 


then we can write 


(4.24) {iedt + - D(X) - U{x))rFGA,oo = 

= ( 2 ^g) 3 d /2 /(^:3;,y,p,g)G^^^p^e‘^"/®G^,p(y)i/>o(y)dydqdp. 

Applying Lemma lD 6 l and adding and subtracting U{Q^), we get 

/ -^n(-^(-iVx + Puf - V{X) - {U{x) - U{QJ) - dtSn)b^Ur,{p,Y) 

+ e{idtb‘^ - ^b‘^d)un{p, Y) 

(4-25) + ed^ ([(Vx 6 " + ib^Pn) + {dtPn - idtQJb^] u„(p, Y)Z-^^ 

^e^bHi [d^Q^Z-^] Un{p.,Y) + [d^ {d^ {b^Un{p,Y)Z-^) Z-'^)] 

+ Pn ■ dtQnb^Un{p, Y) - U{QnWUn{p, Y) 
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Use the Taylor expansion of U{x) about 


(4.26) {U{x) - C7(QJ) = yU{QJ{x - QJ + -y^U{QJ{x - Q^f 

+ ^y^U{Q^){x - QJ3 + lv4[/(Q„)(a; - + E 


|a|=5 


with 

(4.27) 


Rcix) = M _ r)\^\-^d^U{Q^ + t{x - QJ)dT. 

From now on, let us denote the remainder term in (14.261) by R(x,q,p). 
Applying Lemma 14.61 again to ()4.25l) together with (14.261) , we obtain 

/ -n (-i(-iVx + Pn? - V{X) - dtSn)b^MP, 

+ Pu ■ dtQ^b^Unip, Y) - U{QJb^Un{p, Y) 

+ e{idtb^ - ^b^d)ur,{p,Y)+ed^ {VU{QJb^Un{p,Y)Z-^) 

+ ed, ([(Vx&" + i6"Pn) + {dtPu - h„(p, Y)Z-^) 


(4.28) 


+ e- tr [5.Q J/ - V^U{QJWun{p, Y)Z-^] 


+ £2- tr [a.(a.((/ - V^UiQJWur^ip, Y)Z-^)Z-^)] 

+ £^1 tr [5.(a.Q„V3c/(QJ6^h„(p, Y)iZ-^f)] 

+ £^1 tr [5.Q„5.(V3c/(QJ5^h„(p, Y)Z-^)Z-^] 

- £^1 tr [{d,QJ^V^U{QJb^u„ip, Y)iZ-^f] + R{x, q,p)b^Un,p{Y). 
Let us define three operators Lq, L", and L 2 acting on <i> = ^{t, x, y, q,p) by 


(4.29) 


(4.30) 


L 


$ 


and 

(4.31) 


($) :=(-i(-iVx + Pnf - V{X) - dtSn) 

+ • dtQ^^ - C/(QJ4> 

:=(i5t4> - ichd) + a. (V[/(QJ4>Z-1) 

+ a, ([(Vx$ + i$Pn) + (atP„ - iatQ„)$] z-^) 
+ ^tr [a.Q„(/-VV(QJ) 4 >Z-i], 


L^i^) tr [a.(a.((/ - VMQn))^Z-^)Z-^)] 

+ |tr [a.(a.Q„V3p(QJ$(Z-i)2)] 

+ ^tr [a.Q„a.(v3p(gj$z-i)z-i] 


-^tr [idM"vMQ„MZ-^r] . 
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We thus arrive at 

(4.32) {iedt + _ t/(X) - Uix))rFGA,o. 

= .. ^. 3 d /2 // / {L^{b^u^iP,Y))+eL-,{b^Ur,{p,Y)) 

+£^L^{h^Ur,{p,Y)) + R{x,q,p)b^Ur,{p,Y)] p{y)il}o{y) dy dqdp. 

Note that by the choice bn,o = an,oUn{Pn, X), the 0(1) term in the integrand on the right hand side 
of (14.321) vanishes as 

(4.33) LQ{an,oit,q,p)u„{p, X)) = an,oit, q,p){-Hp^ + £;„(P„))m„(P„, X) = 0 
for any a„,o- 

4.3.1. Leading order term bn,o- To determine an,o, we set the order 0(e) term on the right hand side of 
(14.32(1 to zero and get 


(4.34) 


LS(b„,iUn(p,Y)) = -L^(b„^ou„(p,Y)). 


We multiply the equation by Un{Pn, X) and integrate over F; this gives 

(4.35) dtan,o = tr (92P„(Vp^P„)Z“^) - ia„,o-4(P„) -yQ^U - ^a„,o tr(52(5„(VQ^P)Z“^). 
Indeed, by integration, we get (index n is suppressed) 

Jun(Pn,X)(-\(-\SJX + Pn)'^ - V(X) - dtSn)blUn(p,Y) dX 

+ ^{iZ„(P„, X){idtbo - h^d)ur,(p, Y) + u„(P„, X)d^ {VUiQJboUniP, y)Z-^) 

+ uUPn, X)d^ {[(yxbo + iboPn) + (dtPn " idtQM nn(p, Y)Z-^) 

+ II„(P„, X)i tr [d,QM - y^U(QJ)boUniP, Y)Z-^] } dX = 0. 

The perpendicular terms in the 6j’s will now drop out and we can symplify this equation to 

(4.37) - (un(Pn,X),dz ([mn(P„,X)Vp„P„ - VxUn(Pn,X) - iM„(P„, X)P„]aoM„(p, X).^"^)) 

+ - ao^(P„) ■ Vq^P - ^ao)M„(p,X) + ^Oo tr(52Q„(/- Vq^P)Z"^)m„(p, X) = 0. 


(4.36) 


Using (12.111) . we observe that 

(4.38) 

(u„(P„,X), [iM„(P„,X)Vp„P„ - VxUn(Pn,X) - m„(P„,X)P„] • dziaoUn(Pn,Y)Z~^)) = 0. 
Hence, we arrive at 

(4.39) aotr((it„(P„,X),52 • [m„(P„, X)Vp„P„ - VxUn(Pn,X) - m„(P„,X)P„]) Z~^) + 


+ ^iStflo - ao^(P„) • ^Q„U - + ^aoti{dzQ„(I - VqJJ)Z^ ^) = 0. 
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To further simplify the equation, observe that 

[un{Pn,^),d^ • [m„(P„,X)Vp„T;„ - Vxu„(P„,X) - m„(P„,X)P^ 

= l{Un{Pn,X),dzUn{Pn,^)) (Vp„P„ “ P n) “ (u„ (P„, X) , 5^ Vx^n (P«, ft)) 
+ i(9zVp„P„ - dzPn) 

= idzPn {UniPn,X),dp^Un{Pn,X)) {Vp^En - P„) 

- dzPn (U„(P„,X), VpVxMn(Pn,X))p + idzPn{^p„En - I) 


(4.40) 


^ hd,Pn{yl^E^-I)- 


Putting this into (14.391) . we have 

(4.41) iiao tr (d^Pnil - Vp^P„)Z“^) + ^iftao - aoA{Pn) ■ '^qJJ - 


+ 2«o tr(9.QJ/ - VlU)Z-^)\ = 0. 


We arrive at (14.351) finally by noting that 

(4.42) 

Z. Z, Zj 

4.3.2. Next order term hn^i- To characterize we set the order 0{e^) term in (14.321) to zero, we have 


ifltr + ^atr [a;,P„Z„^] =^atr = ^a. 


(4.43) ^ Un{Pn, X) (T”(6„,2h„(p, Y)) + T(((6„,iil„(p, X)) + L”(&„,om„(p, X))) dX = 0. 

Let us first derive the equation for oi. We start with (14.431) written in expanded form 

(4.44) 

^ u„(P„,X)|l tr [d.idMl - V2p(QJ)6o«n(p, Y)Z-^)Z-^)] 

+ |tr [9.(5.g„V3p(QJ6oh„(p,F)(Z-i)2)] 

+ ^tr [d,Q^dA^^U{Q„)boU^ip,Y)Z-^)Z-^] - | tr [{dM^V^U{QJboUn{p,Y){Z-y 


1 , 


+ {idtbi - -bid)un{p, Y) + d, (VP(g„)6iu„(p, Y)Z-^) 

+ {[{Vxbi + i&iP„) + {dtPn - idtQM Mp, Y)Z-^) 

+ ^tr [d,Q^{I -V^U{QJ)b,Unip,Y)Z-^] + (-Pp„ + P(P„)) 62 U„(p, X)} dX = 0. 
Making use of the Hamiltonian flow (13.2|) and the identity ()4.38[) . we arrive at 
-tT^{up^,dz ■ [u(iVP„(P„)) - Vxu - iuPn](ai)) 

+ (idtai - aiA{Pn) -^qU- ^ai) + ao^ tr(a,(a,[(/ - V^^P)Z-i]Z-i)) 

+ aiitr(9.g„(/- V^^P)Z-i) + j^aoir(d.{d.QA^%U{Z-^f)) 


(4.45) 
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Then using (I4.40p and (14.421) . upon simplification we obtain the equation for an,i 

dtOnj = - ian^iAiPn) ’ Vq^C/ + tr (dzPn(Vj,^En)^n^) - tr (^dzQ„(VQ^U)Z~^'j 

( 4 - 46 ) + ^a„.otr ( 5 .(a.[(/- V^^C/)Z-i]Z-i)) + |a„,otr (a.( 5 .g„V^^C/(Z-i) 2 )) 

+ ^an.otr - |a„,otr ■ 

Define the operator Q = Id —n„ where n„ is the projection operator onto the nth Bloch wave, b:^ i 
satisfies Wnb^^i = {un{Pn,x),b^ i) = 0 , and is hence determined by applying Q to Lq Y)) = 

—L^{bn,oUn{p,y))- We obtain 

(4.47) bi.ur^ip, Y) = - (L^)-' Q (&„,o«n(p, Y))). 

Note that the inverse of the operator Lq can be defined on its range. 

Thus, we have obtained the equations for a„^o (14.351) . an,i (14.461) . and b:^ ^ (I4.47|) . This can be continued 
to higher orders. Let us summarize the estimate of these terms in the following propositions. 

Proposition 4.7. For each fc S N, the amplitudes a„p and a„_i, given by (14.351) and (14.461) satisfy 

(4.48) sup Mk[an,o] < oo, and sup Mk[an,i\ < oo. 

te[o.T] ’ te[o.T] 

Proof. By (14.2L (14.31) and (14.41) . we see that the right hand side of (14.351) and (I4.46P are bounded by 
some constants independent of q and p times a„^o and Onp, respectively. An application of Gronwall’s 
inequality yields the result. □ 


Proposition 4.8. For each k gN we have that 


(4.49) 


sup Mk[b^ iUniPn,Y)] < oo. 
te[o.T] 


Proof. The equation for b:^ ^ is given by equation (14.471) . We thus obtain a bound by using the spectrum 

{Um{Pm, •)) 4?(-, Y, q,p)) Um{Pm^ X) 


of Lq. We can write 


(4.50) 


(i^)-'(4’)= E 

m^n 


En{Pn) - Em{Pm) 


Let g = min (lA^d) - \En{^) - E „+i(^)|}. Then for each fc G N, we obtain 

Ce[-7i-.7r]‘i 


(4.51) 


Mk[bi^Un{Pu,Y)]<Mk 


= Mk 


^ ^ (llm {Ptu: ')■! bnp (t, *, p')Un (p, ^)) '^niPni 


m^n 


Unjp, Y) 
9 


an,o{t,q,p){ m ; ') 5 n j ■))l 2 (p) Un(Pn, X) 


m^n 


Hence, by Propositions 14.41 and 14.71 it suffices to control 
(4.52) Mfc ^ (u™(P„,A:),u„(P„,X))p.(p) 

Since /p !««($, x)\'^dx = 1, Bessel’s inequality implies that the above is finite. 


□ 
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4.4. Proof of Theorem 13.11 We will need the following estimate, which is proved in nni Lemma 2.8]. 


Lemma 4.9. Suppose H{e) is a family of self-adjoint operators for e > 0. Suppose ^f(t,e) belongs to the 
domain of H{e), is continuously differentiable in t and approximately solves the Sehrodinger equation in 
the sense that 

(4.53) i£^(t,e) = H{e)'ip{t,e) + 
where C,{t^e) satisfies 

(4.54) \\({t,e)\\< fi{t,e). 


Then, 

(4.55) 


e-‘‘^(^)/"i/>(0,e) -i/^(t,e) 


< e ^ / pl{s, e)ds. 

Jo 


Moreover, for the Fourier integral operator, we have 


Lemma 4.10. If, for fixed x,y € u{x,y,q,p) G L°°(n;C), for each n G N and any t, T^{u) can he 
extended to a linear bounded operator on L^(R'^,C), and we have 

(4-56) l|2^(''^)II.Sf(L2(R<j.c)) < ||'«||ioo(R2d.c). 

Proof. The proof of lemma 14.101 is essentially the same as Proposition 3.7 in and thus is omitted 
here. □ 


We are now ready to prove Theorem 13.II 


Proof of Theorem, \S.1\ Computing ie-^ + -e^V^ — V{X) — U{x) applied to (b^^{t, —, q,p)u„(p, , 

we obtain 


(4.57) 


Y)) = T 



+ (<,2) • 


The expressions for Vnp, Vn,i, and Vn ,2 follows from (14.321) by expanding ¥ and the linearity of Lq, Uf, 
and Ltf. By equations (14.331) and (14.341) . Vn,o and Vn,i vanish. The remaining term 

(4.58) <2 = L'tfffF/uffp,Y)) + R{x,q,p)bll^Un{p,Y). 


satisfies Mfc[u^2] < 00 by Propositions 14.3114.71 and 14.81 Finally, applying Lemma [ 4 . 101 and Lemma 
we obtain the inequality in Theorem 13.II □ 
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